Numerical simulation of water flow around a 

rigid flshing net 



Roger Lewandowski ^, Geraldine Pichot 

'^IRMAR, Campus Beaulieu, Universite de Rennes I, 35000 RENNES, France 
^IFREMER, Technopole Brest Iroise, 29280 PLOUZANE , France 



Abstract 

This paper is devoted to the simulation of the flow around and inside a rigid axisym- 
metric net. We describe first how experimental data have been obtained. We show 
in detail the modelization. The model is based on a Reynolds Averaged Navier- 
Stokes turbulence model penalized by a term based on the Brinkman law. At the 
out-boundary of the computational box, we have used a "ghost" boundary con- 
dition. We show that the corresponding variational problem has a solution. Then 
the numerical scheme is given and the paper finishes with numerical simulations 
compared with the experimental data. 
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1 Introduction 



Recent experimental works [3| show that there are less and less fish in the ocean 
because of intensive industrial fishing. Improvement of the selectivity of fishing nets is 
a major challenge to preserve fishing resources. There are still too many juvenile fish 
and fish with no market value are thrown overboard, leading to a real deterioration 
of the marine ecosystem. Solutions must be found to allow those fish to escape from 
the net when caught. 

Selectivity involves a better understanding of the coupling process between the net, 
the surrounding flow and the fish. Measurements at sea could give some information 
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but they are costly, difficult to perform and not easily reproducible (moving net, 
weather conditions, etc). Therefore, one needs to develop a numerical tool to simulate 
this complex mechanical system. 



The mechanical system made of the elastic net alone in a given laminar uniform flow 
with very simple interaction laws has been studied already, see for instance in 27|, 
and 2J]. A first approach of simulations of the flow around an axisymmetric 



To this point, to our knowledge, no 



n 

rigid net has already been performed in [28 
model exists for dealing with the complex question due to the fish. Finally, there is 
also no work concerning the coupling of an elastic net with the flow. Furthermore, 
it seems that today the numerical simulation of the complete system net/fiow/fish 
does not exist. 



In this paper we deal with the study of the flow around and inside a rigid net in 



the axisymmetric case. Indeed, the code written in [28j cannot be extended to the 
fully 3D case. Therefore, the coupling of the deformation of an extensible net with 
the fluid cannot be considered using this code. Then we have sought a mathematical 
model that we have tested in the axisymmetric case and that can be extended to 
the fully 3D case. We have written the corresponding numerical code and performed 
several simulations to fit the physical constants. Recent investigations have proved 
already that 3D extension is possible and is currently under progress (see 2^). This 
allows to believe that it will be possible in the future to couple our fluid code to an 
elastic code for the net to simulate the system fluid/net. 



Our study starts from experiments performed at the IFREMER's tank of Boulogne- 
sur-Mer (France). A net model rigidified by a resin (see Fig. 1 below) was built and 
velocity components were measured during two experimental campaigns. The first 
one (see [loj) used a Laser Doppler Velocimeter (LDV) technique to get velocity 
components along different profiles. The second one conducted by the second author 
of the present paper made use of a Particle Image Velocimeter technique (PIV) . This 
last campaign emphasizes the locations of turbulent structures in the surrounding 
of the net thanks to instant pictures of the flow. It also gives a good overview of the 
mean flow by averages of pictures. Concerning the velocity proffies, similar shape 
were obtained with the two techniques, except slightly lower value with the PIV. In 
term of accuracy, the LDV technique is much better, that is why the LDV profiles 
were chosen as the reference experimental data to validate our code, for example 
see Fig. 15 to 17 at the end of the paper. It is striking how well the experimental 
velocity data fit with the numerical velocity profiles given by the code. 



The experiments show that the flow we have to simulate is turbulent. Therefore, one 
needs a turbulent model. Yet, we have done simulations by using only the Navier- 
Stokes equations and we did not obtain accurate results. Therefore, we cannot bypass 
the Turbulent model. We have adapted to the present case a classical RANS one 



order turbulent closure model (see for instance [15l, 17 , llal). It is made of an 



equation for the turbulent kinetic energy (TKE) and eddy viscosities functions of 



the TKE into the Navier-Stokes averaged equations. The mixing length has been 
chosen equal to the local mesh size. 



Another important feature of the considered system is that the net behaves like a 



we 



porous membrane. Taking our inspiration in [ll] combined to [2j, and [21 
have modeled the net as a porous membrane by penalizing the averaged Navier- 
Stokes equation with an additional linear term like in the so-called Brinkman Law. 
One considers the net as a fictitious domain and one solves the fluid equations in 
the flow domain as well as in the net domain. However it is an open problem to 
validate mathematically this part of the modelization by using the homogeneization 
theory. We only notice that after a right choice of the permeabilty function K (see 
subsection 3.3) the model yields numerical simulations which fit very well with the 
experimental data. 



The other last important feature of our mathematical model is the boundary condi- 
tions at the border of the computational box. On the lateral boundaries, one impose 
the classical no slip condition. At the incoming boundary, the flow is a given flow. 
The problem is what to do at the outcoming boundary. The natural and classical 
boundary condition should be cr.n = where cr is the strain rate tensor. But as 
observed in [3], one risks artificial eddy reflexions. Moreover, with such a boundary 
condition we are not able to obtain a priori estimates. To overcome this difficulty, we 
have adapted the ideas of [7| to the turbulent case. To do this, we have replaced the 
natural condition by a so-called "ghost condition", the technical condition (8) below. 
This condition becomes the natural one when the flow is laminar at the incoming and 
outcoming boundary (see Remark 4.1). Therefore when observing that far from the 
net the flow remains laminar, we can still take cr.n = at the outcoming boundary. 
This is what we did in the numerical simulations. But we stress that the complicated 
condition (8) is inescapable when dealing with the general mathematical problem. 

Our model is given by the system [(31), (39)] and the assumptions are summer- 
ized by [(23), (30)]. For the sake of simplicity, we have chosen to study the general 
mathematical problem in the 2D case thankfully the axisymmetric case can be eas- 
ily derived, but technical modifications are necessary (see for instance in |9|). The 
existence result stated in Theorem 5.1 is our main theoretical contribution in this 
paper. Uniqueness is an open problem, as well as the general 3D case. 



The numerical scheme uses the finite element method in space, an implicit scheme 
in time for the velocity equation and a semi-implicit scheme for the equation satified 
by the TKE. The parameters settings are defined in section 6.4. As shown at the 
end of the paper, the numerical results fit remarkably with the experimental datas. 

The paper is organized as follows. We start by giving some indications on the experi- 
mental framework, then the modelization is described followed by the mathematical 
analysis. The last part of the paper is devoted to the numerical simulations and the 



numerical results. 



2 Experimental framework 

Experiments have been carried out at the IFREMER center of Boulogne- sur-Mer. 
Velocity profiles have been measured inside and around a rigid resin made model 
built by the Boulogne- Sur-Mer IFREMER team (Fig. 1). This model is like an 
axisymmetric rigid 1/6 scaled cod-end net with diamond-shaped meshes. The end 
of the net is filled with a resin mass modelling a one ton catch of fish and trawled 
with a speed of 1.25 m/s. The net profile as well as the catch geometry have been 
derived from an image processing technique. 




Fig. 1. Model of cod-end net built at IFREMER - Boulogne-sur-Mer 

Note that working on a rigid axisymmetric structure excludes accounting for the 
hydrodynamical forces exterted on the net. Moreover, it restricts the study to an 
axisymmetric geometry. But, at least measurements are possible and mathematical 
fiow models can be tested. 

The model is 1 m long and has an outer maximal diameter of 0.45 m. It is main- 
tained with a frame and set at the bottom of the IFREMER tank. This tank enables 
performance of fiow measurements with velocities between 0.2 and 2 m/s. The esti- 
mation of the velocity to apply in the tank comes from a Froude similitude yielding 
an entrance velocity in the tank equal to 0.51 m/s. 

Hydrodynamical measurements have been performed along several profiles (see Fig. 
2). 

One defines a cartesian reference in the tank, the origin being set at the entrance of 
the net. 

A Laser Doppler Velocimeter (LDV) technique was used to to collect the z and 
y components of the mean velocity (measures are time averaged). The z velocity 
component is the main one we study since it has the direction of the entrance fiow, 
and thus the higher values (see Fig. 3). 



z = -; 
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Fig. 2. Profiles considered of the LDV measures 
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Fig. 3. LDV profiles for the z component of the velocity. 
3 Modelization 

Our model relies on three features: 

(1) Seeing the net, in the fluid point of view, as a porous membrane. The goal is 
then to deflne in which manner the fluid is authorized to flow through the net; 

(2) Directly taking the net and the catch into account in the averaged Navier-Stokes 
equations, which leads to averaged Navier-Stokes/Brinkman equations. This 
way, the boundary conditions at the frontiers of the obstacles are implicitely 
imposed; 

(3) Adding a one equation turbulence model to close the system. 

Our study deals with the mean flow. One can make the assumption that the mean 
flow around the net is axisjTumetric. 
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3.1 Axisymmetric hypothesis 



Assume the cod-end net is embedded in a cylinder full of water. Let us consider an 
axisymmetric deformation of the net (See Fig. 4). As the net is modeled by a porous 
membrane the problem reduces to a 2D one, provided an axisymmetric hypothesis of 
the flow. We admit this hypothesis is a strong one but reasonable in the case of the 
study of the mean flow, since turbulent structures are smoothed by the averaging. 

B 




Fig. 4. Geometry and notations 

In the following, one notes 

• Qw the domain occupied by the water, 

• Gn the net domain, 

• Gf the the fish domain, 

• Gc the domain formed by the frame at the entrance of the net model, 

n = Q^UG, G = GnUGfUGc. 

Using the assumption of an axisymmetric flow and the model of an axisymmetric 
equivalent membrane to describe the net, cylindrical coordinates (O, r, z, 9) are used 
in the simulations. At a fixed value of 9, the mathematical problem reduces to a 2D 
one. The artificial cylinder reduces to a rectangle in the reference (O, r, z) and the 
sides of this rectangle are called Fj, F/ and Fq (see Fig. 4). 



3.2 A membrane model for the net 



Finite elements and finite volumes methods are known to be the common numerical 
methods to compute fluid dynamics. A mesh is built to discretize the fluid domain. 
The difficulty of the netting is that it is composed of a great number of meshes. 



Generating a body-fitted fluid mesh, that is a mesh lying on the nodes and the 
twines of the net, would be far too complex and computer time consuming. Then, 
an exact description of the net would be too demanding in computer resources to 
be conceivable. Another model has to be found. 



In the literature, one finds a model of an axisymmetric membrane to deal with an 



axisymmetric porous structure immersed in a fluid (see [28| ) . 



In [28|, the equations are set on the structure location to express a mass transfer 
in the normal direction to the structure and slip effects in the tangential direction. 
Then, the tangential velocity, denoted Ut-, is set to be governed by Shaffman's law 
and the normal velocity, denoted m„, by Darcy's law. 

This leads to express the velocity components at the wall of the axisymmetric struc- 
ture by: 



Ut = B — , 
on 

(1) 

Un = -KVp, 

V 

where n is the outer normal of the structure, p the fluid pressure, K a permeability 
tensor found experimentally, and B a coefficient dependent on the tangential velocity 
and then deduced from numerical experiences. 

To solve the problem, one builts a cartesian mesh from the geometry of the mem- 
brane, using curvilinear coordinates. The velocity and pressure unknowns are com- 
puted using a finite differences method. 

A drawback of this method is that it is based on a cartesian mesh which is not 
convenient to build and to refine locally in the case of a complex net profile. This 
work then cannot easily be generalized to the case of a 3D deformation of the net. 
One has to find a flow model that allows a future coupling with a moving net. 

Let us keep the idea of seeing the net as a porous medium, as this assumption has 
the advantage of making the numerical programming simpler insofar as twines and 
nodes are no longer taken into account. Then, consider the net and the catch as 
domains with a given permeability. 

As shown in Fig. 4, the domain Gn delimiting the net has a thickness much larger 
than the diameter of net twines (which is typically 3 mm). This idea actually came 
from the analysis of the velocity profiles in the z direction obtained by the LDV 
measurements (see Fig. 3). 



One notices on the LDV profiles (see Fig. 3) inner minima of tlie z velocity com- 
ponent. See Fig. 15-16-17 in the following section for a zoom of each profile. Those 
minima have been noted down (see Table 1). The inner profile of the membrane has 
been drawn thanks to those values. The outer profile is in agreement with the profile 
of the model. 
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Table 1 



This approach avoids a costly mesh generation. However, it comes with the difficulty 
of determining which permeability to apply in the different media. The next part is 
devoted to explain how those media are taken into account directly in the equations. 



3.3 A penalization technique 



The second feature of our model relies on a penalization method that allows us to 



take the presence of the obstacles into account directly in the fluid equations |14| . 
21|. [2I. The method consists in solving "fluid" equations in the entire domain, even 
in the net and catch domains. The net domain is seen as a porous medium, and 
the catch domain as a solid medium, where a no-slip boundary condition should 
hold. Those media are explicitly included in the fluid equations by the addition of a 

penality term of the velocity, namely ^ u. This leads to Navier-Stokes/Brinkman 

equations. Notice that such laws have been derived from an homogeneization process 
in other situations, as in [l|. This theoretical question remains open in this particular 
context. The function i^(x) varies from one domain to another. It is a permeability 
parameter that is very small in the solid domains, e.g. the catch, to force the velocity 
to be zero, and very high in the fluid domain, so that the averaged Navier-Stokes 
equations hold and are set to a defined value or function in the porous domain (here 
in the net domain) depending on its permeability. 



At a first glance, the function -ft'(x) is set to be constant by parts. The net domain 
is decomposed in three parts, G^, i=l, 2, 3 (see Fig. 5) of constant permeability 
that is all the more important as we are closer to the catch (see Part 6). In a future 
work, we will try to make it depend on the mesh opening, the mesh angle between 
the mesh and the local flow. 




Fig. 5. Decomposition of the net domain - Notations 
3.4 Addition of a turbulence model 

The third point comes with the average of the Navier-Stokes/Brinkman equation, 
since Direct Numerical Simulation would not be able to treat a problem with such 
a high Reynolds number (here Re = 10^, using as reference length the maximal 
diameter of the catch, i.e. 0.45 m, and the entrance velocity as reference velocity 
that is equal to 0.51 m/s). A kind of Reynolds Averaged Navier-Stokes (RANS) 
turbulence model is then added to close the system of equations. It consists of one 
equation for the turbulent kinetic energy. The averaged NS/Brinkman equation and 
the turbulent kinetic energy equation are coupled by the means of a eddy viscosity, 
denoted Uf 



4 Description of the mathematical problem 

4-1 The domain 

We return back to the description of the geometry. As already said, the flow under 
study is axisjTiimetric. In order to avoid technical complications, we have chosen to 
study the mathematical problem set in a domain in IR^. We refer to [sl to go in 
further developments in the axisymmetric case. 

The boundary F of the computational box is defined by the input board Fj, the 
lateral boards F; and the artificial output board To, 

Ti = [0,A], = (0,0), A={0,a), 

Ti = [C, O] U[A,B],B = iP, a),C= {(3, 0), 

ro = [B,c], 




f> If 

Fig. 6. Description of the domain 



4-2 The equations 



The unknowns are : 

• the mean velocity vector field u = u(t,x) = [u^ [t , x.) , [t , x.)) , x = 

• the mean pressure scalar field p = p{t, x), 

• the turbulent kinetic energy k = /c(t,x) 



'x,y), 



One defines the deformation tensor s by 



£(U) 



Vu + Vu^ 



(3) 



The turbulent strain rate stress tensor a is defined by 



<t(u,p, k) — 2 Vtik, x)£(u) — pid. 



(4) 



The Reynolds Averaged Navier-Stokes turbulent closure model of order one including 
the Brinkman laws, is given in [0, T] x (T > 0) by the following equations, where 
£ > is fixed. 



dtVi + (uV)u - V • o-(u,p, k) + 0(Ig,uGc) + J^^Gn + £ln. j u = 0, 
V • u = 0, 

dtk + u.Vk-V ■ {nt{k,x)Wk) = 2iyt{k,x)\s{u)\^ - S{k,x}. 



(5) 



In the equations above, and are the eddy viscosities and £ the backward term. 
Their analytical expressions are given in section 4.4 below. 



4-3 The boundary conditions and the initial data 



4-3.1 Boundary conditions 

The input field Ui = (mj, 0) at the boundary Fj is a data of our problem. The 
boundary conditions we consider are the following. 



on Fj : u = Ui = (mj, 0), k = 0, (6) 
on F; : u = 0, k = 0, (7) 

<t(u,p, A;).n = -^(u.n)"(u - Ui) + (u.n)ui 
on Fo : <( ^ (8) 

k = 0. 



In the formulae above, Ui denotes the field equal to {ui{x — /3,0)) on Fq. One uses 
the boundary condition (8) in order to avoid eddy reflections at the open boundary 
Fo and to be able to prove the existence of a dissipative solution to the system (5). 

Remark 4.1 The natural boundary condition for the velocity at Fq should be 
(t(u,p, A;).n = 0. In f^], the authors study the case of the Navier-Stokes equations 
without a turbulence model and in a channel without a rigid body. They remark 
that the boundary condition a.n = yields numerical eddy reflexions at the out 
open boundary. Moreover, the existence of a dissipative solution is not known in 
such a case because of a term /py(u.n)|up which appears in the energy equality due 
to the convection. Without additional information on the sign of (u.n) at Tq, no 
a priori estimate is avaible. This is why the authors in change the boundary 
conditions. We also change the boundary conditions by an adpatation to the case 
of our turbulence model. Notice that when the flow is laminar at Fq and (u.n) > 
on Fq, the boundary condition reduces to the classical one up to the term (u.n)uj. 
This is an additional forcing term. Without this term, it is easy checked that one 
can only derive an a priori estimate when a smallness assumption on Uj is satisfied, 
an assumption which would restrict the problem to a laminar one. Therefore, this 
term seems to look coherent when the flow is turbulent at the incoming boundary. 
However, in the numerical simulations we have taken a{u,p, k).n = 0. Indeed, the 
experiments suggest that the flow is laminar far from the net. Therefore, our choices 
are in concordance with reality and yields a rigorous mathematical analysis. 

Remark 4.2 For convenience and the sake of simplicity, we have chosen to develop 
the theoretical part by fixing k = at T^. A more natural boundary condition at 
Fo zs = 0. This is the condition that we use in the numerical simulations. 

From the mathematical viewpoint, we then have to impose = — (u.n)^/;; at 

To. Therefore the discussion in remark 4-1 above holds in this case. However, this 
boundary condition yields serious mathematical complications that would have been 
out of the scope of this paper. In subsection 5.6 we give some explanations about this 
case. 



Throughout the paper, we assume that 



(9) 



4.3.2 Initial data 

The initial data are specified by 



Vx e n, u(0, x) = Uo(x) G iL\Q)y, (10) 

Vxel], A;(0,x) = A;o(x) e Li(l]). (11) 

Moreover, we shall assume that Uq satisfies the compatibility conditions 

V • uo = 0, (12) 

UQ.n = Ui on Fj, (13) 

Uo.n = on T^. (14) 



Remark 4.3 The assumption uq G (L^(r2))^ gives a sense to u^.n in the spaces 
(ifQO^(rj))' and {HQQ^(Ti)y, making (13) and (14) consistant as a consequence of 
(9). 



4.4 The eddy viscosities and main terms 



4.4-1 Eddy viscosities 

The eddy viscosity function Vt is a non negative bounded function of k and x 
equal to uq + £(x)y^r + when G [0, /cj for a given fee and p > is fixed. The 
viscosity vt is thus given by 



i/t(/c,x) = i/o + £(x)Y'r+ when |A;| < /cc, (15) 
i't{k,:x.) — V2, when \k\ > kc+1, (16) 

o N f K^) . . \ ,3 //(a;)(-3A;,-2) 



2Vt + /cc 

kl - 3kl) when /c^ < ^ < + 1 

(17) 



2 V'?' + kc 



where r > 0, kc> 0, V2 > Vi — 1/0 + £{y:)\^T + kc 



Fig. 7. Shape of vt 

The function £(x) is a local scale of the flow. It is a non negative bounded function 
of X on with 

Vx G f], < £o < <Lq<oo. (18) 
The eddy diffusivity /if is of the same form as Vt and 



z/Q + C£(x)^r + 1^1 on the range [0, (19) 
for C > and f > flxed coefficients. 

4.. 4- 2 Backward term 

The backward term S{k,x.) is given by the formula 

S{k,^) = -^k^/k. (20) 

4 -4 -3 Permeability 

The permeability function K[x) is a continuous function that satisfies 

yxen, 0<Ko< K{x) <Ki<oo. (21) 
In the remainder, one shall set 

P(u)(t,x) = (^^(]Ia,uG.(x)) + Y(^^gA^) + eln^^ u(t,x), (22) 

where e > is fixed. 



5 Mathematical analysis 



5.1 Main result 



We summerize the hypotheses: 



vt e C\ V (A;, x) e IR X fi, Q<vq< i^tik, x) < N < oo, (23) 

Ht&C\ V (A;, x) e IR X < Ho < T^t{k, x) < M < oo, (24) 

^ e Vx e < 4 < £(x) < Lo < oo, (25) 

S{k,x) = j^k^l (26) 

G C\ Vx G 17, < i^o < i^(x) < i^i < oo. (27) 

uo G L^(n), V-uo = 0, Uo.n|ri = Mi, Uo-nlr^ = 0, (28) 

A;oGL1(Q), A;o>Oa.e, (29) 

u, e i/oo^'(r,). (30) 
The problem is the following 

dtu + (uV)u - V • (t{u,p, k) + P(u) = 0, (31) 

V • u = 0, (32) 

dtk + n .Vk -V ■ {nt{k,yL)Vk) ^ 2vt{k,x)\e{vi)\^ - £{k,x). (33) 

Vxel^, u(0,x) = uo(x), (34) 

VxeQ, fc(0,x) = A;o(x), (35) 

u|r, = Ui = (mi,0), A;|ri=0, (36) 

u|r, = 0, k\T, = 0, (37) 

<r(u,p. A;). n|r„ = --(u.n)-(u - Ui) + (u.n)ui, (38) 

fc|r„ = 0. (39) 



Our main result is the following. 

Theorem 5.1 Assume that hypotheses [(23). ..(30)] hold. Then Problem [(31). ..(39)] 
admits a solution (u,p, k) on any time interval [0, T] in the sense of the distributions, 
where 



u G L2([0, T], {H\n)f) n L~([0, T], L2(0)), (40) 
PGL2([0,T]x1]), (41) 

k e L^{[Q,TlL\n))n{{^ L^{[Q,TlW''^{m- (42) 



Remark 5.1 Uniqueness remains an open problem. 



5.2 Lifting the boundary condition 
5.2.1 Auxiliary Stokes Problem 

In this section, we describe how to hft the boundary conditions to reduce the problem 
to a problem with homogeneous boundary conditions on FjUr;, as it is usually done 
in mathematical problems where Navier-Stokes Equations are involved. 

Recall that Vt^, is the water domain and G the net domain (see section 3.1). 

The incoming flow Uj is prescribed at the boundary Fj. We define Uj on the output 
boundary To and still denote it by Uj, the field defined by 

Vx = {x,y) G To, Ui(x,2/) = {ui{x - /9,0)). 

Let us consider the Stokes problem 

-Avo + Vgo = in 

V-vo = inl]^, (43) 
vo = g on r U dG, 

where VL^ is the water domain, G the domain delimited by the net and g is the field 
defined by 

on Vi U To, g = Ui, 
on Ti U 9G, g = 0, 
Notice that the following compatibility condition is satisfied: 

j g.n = 0. (45) 

In the following, we note 

Ll{n) = {qeL'^{n)- j g(x)c;x = 0}. 

Theorem 5.2 Assume thatUj G HQ^'iVi) (assumption (30)j. Then Problem [(43) — 

(44) ] has a unique solution (vo,go) ^ H'^iS^w) x {H^iS^w) H Lg(fi^)). 

Proof. On one hand, it is established by Corollary 5.9 in [sj that g G [H^^'^iVyf' 
because Ui G HQ^'iVi). On the other hand, g satisfies the compatibility condition 

(45) . Moreover, fi^, is a convex polygon in dimension 2. Therefore, applying Theorem 
5.4 and Remark 5.6 in jll| §5 (see also [l2|)) one knows the existence of a unique 
(vo,go) e H'^{^w) X {H\VlJ) n Ll{VL^)) solution to Problem [(43) - (44)]. 



(44) 



Remark 5.2 In practical situations, Uj is a Poiseuille flow. Therefore, one has 

Uj{x, y) = Uj{y) = Dy{a - y), 
where D is a constant. We first note that U[ is C°° on Fj. Moreover, one clearly has 

/ ay < +00 ana / -ay < +00. 

Jo y Jo [a-y) 

Therefore, thanks to the definition of H^q'^ (see in fs^], chapter 1, §ii or in [8] chap- 
ter 6), Ui & ifoo^(rj). Unfortunaly, Ui ^ ifoo^(rj). Therefore, one cannot guaranty 
that g G iJ3/2(r) 

using the results above mentioned and only g G H^^'^iV). In 
such a case, only regularity for the velocity can he obtained a priori and that is 
not enough regularity for what follows, as we shall see in the remainder. 

Remark 5.3 Since g G [i/^/^(r)]^, the trace on To of eiyo) is in [H^/'^{Vo)Y well 
as the trace of qQ on To is in H^^'^iVo)- Then, because Vt is a bounded function, for 
every k G L^{Q), 

a{^ro,qo,k)e[H'/'iTo)]\ (46) 

From now, one still denotes by vq the field defined on whole Q and equal to vq in 
flu, , the velocity part in the solution to Problem [(43) — (44)], and equal to inside 
G. Since 

. H\n^) c L-(n^) 

• dG is of class C^, therefore one can use Proposition IX. 18 in [g], 



one has 

voe H\n)nL°°{n) (47) 

and 

where C only depends on a and /5. By extending go by zero outside fl^ and still 
denoting the expension by q^, one has 

CTivo,qo,k)e[L'in)]\ (49) 

Notice also that 

P(vo) = 0. (50) 

5.2.2 Change of variable 
We set: 

u = u + Vo, p = p + qo. (51) 



It is straightforward to prove that (u, p, k) is governed by the following system: 



atU+(uV)u- V-o-(u,p,A;) +P(u)+ 

< 

^ (uV)vo + (voV)u + (voV)vo - V • <r(vo, ?o, ^) = 0, 
V • u = 0, 

dtk + u .Vfc - V • (/it(A;, x) VA;) = 2vt{k, x) |£(Ci) 1^ - E{k, x) + 
4z/i(/c,x)£(Ci).£(vo) + 2z/i(A;,x)|£(vo)|2 - voVA;. 

u|t=o = uo - vo, k\t=Q = ko, 
u|r,ur, = 0, klviUTt = 0, 

cr(u,p,/c).n|r„ = 

-2[(^ + Vo)-n]~u + [(u + Vo).n]vo - <r(vo, ?o, ^)|ro-n, 
A:|r„ = 0. 

5.5 Variational formulation 
5.3.1 Functions space 

The natural space for studying Problem [(52) (57)] is the space 

V^{we{H\n)f- V-v = 0; v|r,ur. =0.} (59) 

In order to use De Rham Theorem and have an Inf-Sup condition on the pressure, 
we must check that smooth vector fields with null divergence and equal to zero on 
Ti U Ti consitutes a dense space in V . This is the goal of what follows. 

Let B = (2/9, a), C = (2/3,0) and let (l be the square in bounded by the 
points O, A, B and C. Let s be the symmetry through the axis x = (3, that is 
s(x,y) = (2(3 -x,y). 

We also denote by the square bounded by the points C, (7, B and B, also defined 
by = s{n). 

Let V be the set 

V = {v e iH\n)f; V-v = 0; vj^^ = o} 

as well as 

V = {v e (P(Q))'; V-v = 0}. 
Being given v e V, let be its restriction to the square Q. One obviously has 

e V. 



(52) 

(53) 

(54) 

(55) 
(56) 

(57) 

(58) 



Being given v ^ V, let v'^ be its extension to Q defined as follows: 

Vx G n, v^(x) = v(x), 
Vx e v^(x) = v(s(x)). 



(60) 



Notice that G and one has 

/ |Vv^|2 = 2 / iVvl^ VpG[l,oo[, / |v"|P = 2/|vP. (61) 
Jfi" Jn Jn" Jn 

Finally let V be the space made of the restrictions to ^2 of fields in V, which means 

V = {v G [C°°(^])]; 3 V G V s.t. V = v^} (62) 

We prove the following. 

Lemma 5.1 The space V is dense in V . 

Proof. Let v G V^. Since Cl is simply connected and has a Lipschitz boundary, one 
knows thanks to Corollary 2.5 in ll| that V is dense in V. Therefore, there exists a 
sequence (w„)„£]n of fields in V that converges to v*^ in the space V. One obviously 
has 



/ |V((w„),-v)|2< / |V( 
Jn Jn" 
This shows that the sequence ((w„)r)ri6]N converges to v in and each (w„)j. lies 
in V by definition. The lemma is proven. 



5.3.2 The variational Problem 

For the sake of the simplicity, up to now and throughout the paper we shall note 
I'tik) instead of ^'t(/c, x). Notice firstly that V (vi, V2) G V^, V (/c, q) G V^fl)"^ one has 

- / (V ■ cr(vi, q, k)).Y2 = - (cr(vi, g, k).n).\2 + / 2 Vt{k)£ (vi) : ^(va). 
Jn JTo Jn 

The variational formulation of the problem is the following, where the pressure does 
not appear anymore and will be recovered using The De Rham Theorem. In the 
following, one denotes 

W = L2([0,T],F) nL°°([0,T], {L'iQ))^). (63) 



Find 



u G VF, u(0,x) = Uo(x) — Vo(x) a.e in Vt 

keL^{[o,T],L\n))ni n L^i[o,Tiw'^^m) 

p<4/3 



(64) 
(65) 



WITH 

dtu G {[0,T],V')nW\ (66) 

AND SUCH THAT V V G L^([0, T], \^), 

<dtU,v>+[ [{iiV)u.v+[ j 2ut{k)e{vi):e{y)+ I [ -[(u + vo).n]-ii .v 

^ Jo Jn ^ Jo Jn ^ Jo JVo ^ 

[ / [(ii + vo).n]vo.v+ / [v{u).v+[ / [(voV)(ii + vo).v+ 

Jo JFo Jo Jn Jo Jn 



r [ 2z/i(fc)£(vo) :£(v)+ r [ (o-(vo,go,A:).n).v = 0, 
Jo Jn Jo JVo 

(67) 

FOR ALL r G P'([0,T] X Vl), WITH r(T, ■) = 0, 



r / 5tr/c+ / ((ii + Vo)V)A;.r + [ fitik)Vk : Vr + [ A;o(x)r(0, x)rf> 
Jn Jo Jn Jo Jn Jn 

[2u,{kMu)\' - £{k) + u,{k){Ae{i,)e{w,) + 2|£(vo)n]r 



T 

Jn 



(68) 



5. 3. 3 Consistency of the variational formulation 



The variational formulation for the fc-equation is the classical one, as in [16|, [17 



18| and [15| . The variational formulation for the velocity is also classical up to the 
boundary terms. Each boundary term where vq is involved is nice since vq does not 
depend uppon the time and is equal to Uj on To which is in particular in L'^{Ti). 



Nevertheless the term 

rT 



[u.nj u.v 

'0 JFo 

is fearsome. We prove the following lemma which guarantees the consistency of the 
variational formulation above. For the sake of simplicity and as far as no confusion 
occurs, we still denote by v the trace of v for any v G W. Moreover, one defines the 
W norm by 

\W\\w = ||v||L2([0,T],y) + ||v||loo{[o,T],(L2{Q)2)- 

Lemma 5.2 Let (u, v) e W x W . Then 

rT r 

u.nl^u .V < C||u||L||v||m/, (69) 

JFo^ 

where C is a constant that only depends on a and (3. Moreover, there also exists a 
constant C such that 

VvGL«/^([0,T],l^), f I [ii.n]-u.v<C||u||2^||v|L8/3([o,T]y) (70) 
Proof. Let v G W . On starts from the classical interpolation inequality (see in |20i] ) 

II II ^11 llV4|| ||3/4 

I I V| 1^3/4 < ||v||/2 ||V||^. 



One deduces that 

I I V| 1^8/3(^^3/4) < ||v||J;'i(^2)||v||^'l'J^) < -||v||Loo(i2) + -||v||L2(y) < \\v\\w (71) 

One deduces by the trace Theorem that 

l|v||L8/3(Hi/4(ro)) ^ C'||v||iy. (72) 
Moreover, thanks to the Sobolev Theorem, 

l|v||L8/3(i4(r„)) < Cllvllw'. (73) 
Let u e ly. It is clear that at Tq, 

u.n e L~(//-V2(ro)) n L''{H'/''{ro)). 

By using again a simple interpolation inequality one deduces easily that 

||u.n||i4(i2(r„)) < Cllullw'. (74) 
Therefore, (u.n)u e L^^^{L^^^{Tq)), as well as (u.n)"u and one has 

1 1 (u-n)"u| |L8/5(L4/3(ro)) < C| |u| 1^. (75) 
The rest of the proof is now a direct consequence of (71), (73) and Holder inequahty. 

5.4 A priori estimate 



Proposition 5.1 There exists a constant Ci = Ci{uQ,Ui,i>,a, f3) and for each p < 
4/3 a constant C2 = C'2(uo, Mj, u, fi,p, a, P) such that for any smooth solution (u, k) 
to the variational problem [(68), (67)] one has 



u||L2([o,T],y) + ||u||loo([o,t],l2(q)) < Ci, (76) 
< C2. (77) 



Proof. We proceed in two steps. We first estimate the velocity and then the Tur- 
bulent Kinetic Energy (TKE). 

Step 1. Estimating the velocity. One multiphes the equation (52) by u and integrates 
on Q. A technical but easy computation using the boundary condition u (57) yields: 



1 ci 

-— u 

2dt " Jn Jn ^ (78) 

I V{n).VL- I Vo(8)(u + Vo) : Vu+- / ((u + vo).n)+|u|2 = 0. 
Jn Jn 2 Jto 



Since 

0<fv(u).u and 0<-/ ((li + vo).n)+|iiP, 

Jn 2 Jr,) 

using (23) and (47), the energy equality (78) yields 

^^l|u||i2(f,) + ^2z/i(A;)|6(u)|2 < iV^|£(ii)||£(vo)| + ||vo||ooL|fi||Vii| 

+ l|vo||L / |Vu| 



(79) 



By using Young and Korn's inequalities, one has 



,^(ii)lk(vo)|<^/ k(vo)P + |/ \e{u)\^ (80) 
n zQ Jn 2 Jn 

|fl||Vii| <\l- + c\ f |u|2 + ic f |s(u)P. (81) 



n \2( J Jn 2 Jn 

where ( will be fixed later on and C is the constant in the Korn inequality. Finally, 
by always using the Young inequality combined with the Cauchy-Schwarz inequahty, 

INlL/jva|<|||v„||i4/j.(a)P. (82) 

Thereofore, (79) combined with (23) yields 



|u| li.(o) + (^0 - (C/2)(iV + C + 1)) \s{n)\' < + llulli. 



+ ^ / k(vo)P (83) 



2C Jn 



af3 
+^l|vo 



|4 
I oo 



We choose C be such that (z/q - (C/2)(A^ + C + 1)) = uq/2. One deduces from (83) 
and Gronwall's lemma, combined again with Korn's inequality, the existence of C = 
C(ui, N, uo, a, P, T, uo), which blows up in a rate and such that 

W^Ww = ||u||L°°([o,r],L2(n)) + ||u||l2([o,t],v) < C. (84) 



Step 2. Estimating the TKE. Notice first that by using the same arguments as in 



16 



or m 



ITj, one can make sure that k > a.e. as far as we assume ko > 0. The 



boundary terms does not create any troubles because 



r [ iiu + Vo).n)-ki-k~)= r [ ((ii + Vo).n)-(fc-)2>0. 
Jo JFo Jo JFo 



The other terms are hke in the general situation studied in [l7j chapter 4. From now 
and throughout the rest of the paper, one works with k > 0. 



Thanks to (84), we can use the Classical Boccardo-Gallouet estimate (see [^). By a 
proof already done in 16j, 17|], 18| and 15|] and since we are working in a 2D case 
and A; = on dQ, one deduces that 



3C = C(ui,iV, z/Q, a,/?,T,uo); | |L°°([o,T],Li(n)) <C, 



i5) 



and 



Vp<4/3, 3C = C{p,u„N,uo,a,p,T,UQ); \\k\\LP^[o,T],w^'P{Q)) ^ (86) 



5. 5 End of the proof of the main Theorem 



The proof now is the same as in 16|], 17|], 18|] and 15|, up to the additional terms 



due to the extra boundary conditions for the velocity. We construct a sequence 
of smooth approximated solution (u„, A;„)„gN (for instance by troncating the l.h.s 
of the fc-equation and using the Galerkin method). The trick is to prove the weak 
convergence in L^{[0, T], V) of the sequence (u„)„g]N (up to a subsequence) to u G 
which satisfies the formulation (67), and in particular, that can be taken as a test 



function in (67). Once this task is finished, the rest is classical and works as in 16 



17] , 18| and [l5| since we already have obtained all the required a priori estimates. 



Let u G W. One has, after a part integration on the convective term, 

<c)tUn,v>- / /u„(g)u„Vv+ / / (u„.n)u„.u+ / / 2ut{kn)e{\ir, 
^ ^ Jo Jn Jo JVo Jo Jq 

+ / -[(u„ + Vo).n]"u„.v - / / [(u„ + Vo).n]vo.v+ 

Jo JTo I Jo JVo 

r I P(ii„).v+ r ( [(voV)(ii„ + vo).v+ r / 2ut{k)e{w^):e{Y) + 
Jo^Jn Jo Jn Jo Jn 

/ / (cr(vo,go,/i;)-n).v = 0, 

Jo JFo 



17) 



By using (84), one knows that the sequence (u„)„gN is bounded in W and one 
may extract a subsequence (still denoted by the same) that weakly converges in 
L^([0, T], V) and in L°°{[0, T], L^) to some u G W. One needs compactness, and for 
it we shall use the Aubin-Lions Lemma. Of course, all the terms involved in (67) 
satisfied by each u„ are nice except the terms 

u.„.n)u„.u and / / [u„.n]"u„.v (88) 
JFo Jo JFo 



which are the worse terms and which constitutes the only new difficulty in this 
problem compared with previous works already quoted. Thanks to inequality (70) 
combined with (84), the applications 



/ / (u„.n)u„.v, 



u„,n u„, .V 



10 



are bounded in the space L^/5([0,T],\/') ((84) holds for the second one, the proof 
is the same for the first one). Since we are working in a 2D case, and thanks to 
the regularity of vq, all the other terms are bounded in L'^{[Q,T],V'). Therefore, 
the sequence {dtUn)n€i^ is bounded in L^/^([0, T], y) as well as in W. Applying 
the Aubin-Lions Lemma, one concludes that the sequence (u„)„g]s[ is compact in 
L*/^([0,T], (L^(fi))^). Hence we are back to the usual situation concerning compact- 
ness in this type of problem. We bypass the details. We still denote by (un)nGN 
the subsequence which converges to u almost everywhere in fl and strongly in 
L^([0,T], {L\n)y) (we are the 2D case). 



One has analogous compactness properties for the sequence {kn)n€i<i which converges 
weakly in each -^^([0, T], Wq'^{^1)) (up to a subsequence andp < 5/4) to some k in the 
space np<5/iLP{[0, T],Wo'^{n)), almost everywh ere in Q and stronly in L'^{[0, T] x Q) 
for some q > 1. 



Passing to the hmit in all the terms in (87) is a classical game and follows proofs 
done already in previous papers (we are in the 2D case) , except concerning the terms 
(88). We show how to pass to the limit in the first one, the second one being treated 
by the same reasoning. Notice that one has 

C H^^^ C V, 

the injections being dense and compact. Hence, the sequence {un)neN is compact 
in the space L^/^([0,T], (i7^/*^(f2))^). By uniqueness of the limit, it converges to u 
in this space. Following the chain rule of the proof of Lemma 5.2, one deduces that 
(u„.n)„g7v converges strongly to u.n in L^([0, T], L^(ro)) while (Un)n6N converges 
strongly to u in L^/^{[0,T], {L\ro)y. Therefore, 



lim / / (u„.n)u„.u= / / (u.n)u.u. 



The rest of the proof is now classical. 



5. 6 Neuman Boundary Condition Type for the TKE 



We are now working in the case where k does satisfy on Fq 



dk 

f^tg^ = -(u.n)-/. (89) 



instead of k = 0. Because this case yields serious mathematical complications, we 
shall not give a complete proof of the existence result. We shall limit ourself to 
locating the difficulties, giving the main a priori estimate and to indicating the 
direction to take. Details will be written in a forthcoming paper. 



5. 6. 1 Variational Formulation 

When k satisfies (89) at Fq instead of /c = 0, the variational formulation for the 
k-equation becomes: for all r G C°°([0,T] x Vt), with r|r,ur, = and r(T, ■) = 0, 



-/ f dtrk+ f /((ii + Vo)V)A:.r+ / f fit{k)Vk : Vr+ 
JO Jn Jo Jn Jo Jn 

[ f {{u + Vo).n)-kr+ f ko{^)r{0,^)d^= (90) 
Jo__Jro Jn 

[2z/,(A;)|5(u)p - S{k) + u,{k){Ae{u)e{^r,) + 2\e{^ro)\')]r 



T 

Jn 



The source of difficulty is the additional term 

4 = r f ((ii + vo).n)-/cr. 

Jo JFo 



5.6.2 A priori estimate 

One starts first with the a priori estimate. We show that in the following, there is a 
situation where the Boccardo-Gallouet result [4] can be applied. 

Let g be any non decreasing non negative piecewise bounded function defined 
on IR"*", G{k) = g{k')dk' . Notice that G is non negative and thanks to the mono- 
tonicity of g, one has 

VA;GlR+, 0<kg{k)-G{k) (91) 
Therefore, by choosing g{k) as test function in (68), with u = u + vq, one has 

^G{k)+ I iimg\k)\yk\^+ f (u.n)+G(fc)+/ {u.^Y {kg{k) - G{k)) = 
at Jn JFo JTo 

f g{k)[2u,{kMu)\' - m] 
Jn 

(92) 

Since g is non negative, combining (47), (84) and (91) one has 

'^-G{k)+ [ fit{k)\g'{k)\Vk\^ <d\\g\U d = d{u„N,iyo,a,(3,T). (93) 



dt Jn 



Therefore one can deduce that the results in [3] apply. Therefore the estimates (85) 
and (86) still hold in this case. 



5. 6. 3 Consistency of the variational formulation 

As said already, the difficulty is due to the term Ik- Recall that from the proof 
of Lemma 5.2, ((u + vo).n)~ G L^([0, T], L^(ro)). On the other hand, by com- 
bining the trace therorem with the Sobolev Theorem, it easily checked that k G 
np<4/3^^([0,T],L^(ro)). Here the critical case is the space L^/^{[^,T],L'^{Tq)), 
which is not achieved. Therefore, it is not guarantied that the integral Ik is defined. 



The way to go round this difficulty is to renormalize the equation for /c, as in [17 
chapter 5 and also in [2^. Roughly speaking, one does not take a test function r in 
the equation, but rip{k) for functions ip having compact support. Then Ik becomes 



4,v = / / ((u + Vo).n) ktp{k)r, 

Jo JFo 



which is defined since kiplk) is bounded. Of course, when doing this, new terms 
appear in the variational formulation. This is now out of the scope of the present 
paper and will be the subject of a next paper. 



6 Numerical simulations 



Simulations have been performed using the free software Freefem-|--|- (see [13|). It 
allows computations of 2D and axisymmetric fluid dynamics by the means of the 
finite elements method (FEM). 



Remember that the net is modeled as a porous membrane and enclosed in a fictive 
cylinder. Assume that flow is also axisymmetric. Recall this is a strong hypothesis 
but reasonable in the case of the study of the mean velocity around a rigid net. 
Then the problem reduces to a 2D one. The geometry shown on Fig. 8 and drawn 
in Freefem++ has an outer net profile and a catch profile in agreement with the 
model of Boulogne-Sur-Mer. The inner net proffie is defined by the minima of the 
z component of the velocity located on the LDV profiles (see table 1). To take into 
account the difference of permeability of the net (mainly due to the variations in the 
mesh opening), the domain Gn has been decomposed in 3 sub-domains: G^, G^, G^. 

Let us work in cylindrical coordinates, the z axis being the revolution axis of the 
membrane: 



X 


= rcos9, 




' y 


= rsinO, 


(94) 


z 


= z. 





m 



Fig. 8. Geometry of the net 

Let n = {{r,Z,9),r e [rmin, rmax], Z G [Zmin, Zmax], d ^ [0, Tt]}. 

Let u = {ur,U0,Uz) denote the mean velocity unknown in cylindrical coordinates. 

Assuming a planar flow, then uq = and thanks to the axisymmetric hypothesis, 
derivatives with respect to the variable 9 are zero. 

At a fixed value of 6, we work on a 2D domain: 

Notice we keep the notations: fi^ for the fluid domain, Gc for the ring that maintains 
the model inside the tank, for the membrane (net) domain and Gf for the catch 
domain. 

In the following, the operators (gradient, divergence,...) are considered in cylindrical 
coordinates. 

The solid part has a very small permeability, denoted Ks{r, z) <^ 1, leading to force 
the velocity to be zero in that part (then forcing a no slip boundary condition). 

The porous part has a permeability chosen here to be constant by subdomains G^, 
denoted Kqi^, i = 1,2,3. The fluid domain has an infinite permeability so that the 

penalization term vanishes in that part, denoted Kf{r, z) ^ 1. 

The coupled problem (31) — (39) is implemented under the following variational 
form. 




6. 1 Weak formulation 



At first, let us assume that there is no reflexion at the outer boundary and consider 
the boundary conditions (38) — (89) reduced to: 



<r(u,p,fc).n|r„ = 0, (95) 
^k = 0. (96) 



Moreover, let us replace the no slip boundary condition for the velocity (see equation 
37) on Ti by slip boundary condition and the homogeneous Dirichlet condition for 
k on by a non homogeneous one: 



= 0, = 0, k^kQ sur T,, (97) 

or 



Denote V{ytr,z) and Q{ytr,z) the space defined as: 

V{Vtr,,) = {v e {H\nr,z)f: V, ^ sur T^, Vr ^ sur Ti U r^}, (98) 



Q{^r,z) = {g e L\Qr,z)}. (99) 



and 



yV{^r,z) = {w e L'^{nr,z),w = sur U Ti}. (100) 



A weak formulation of the coupled problem (31) — (36) , with the boundary conditions 
(95) - (96) - (97) yields: 



Find (u = {ur-, Uz),p, k) e V{Q:r,z) x Q{^r,z) X yV{Q:r,z) such that: 



I —\'\r\Tidrdz-\- j (uV)uv |r|7r(ir(i2; — j pV ■ yr\r\'ndrdz 
+ i / (i/Q + i^t)(Vu + (Vu)*) : (Vv + (Vv)*) \r\'Kdrdz 
+ V{u){t,{r,z))v\r\'Kdrdz 

(101) 

- /" V ■ uq\r\T:drdz = eV{VLr,,),\lq e Q(fir-,^); 
/ —w\r\'Kdrdz-\- j {\xV)kw\r\'Kdrdz + ( Vti^k •.Vw)\r\'Kdrdz 

— / ^|Vu + (Vu)*|^i(; |r|7rcirci2; + j k^w \r\-Kdrdz — 0, 

with 



(102) 



- (^(Ic,uc.(r,.))+g^-^I.,(r,.) 



^.^ Finite elements discretization 



Using the mesh generator of Freefem++, one builds an unstructured mesh 7^ of the 

domain -[(r, ^),r G [Tmim'^max\i ^ ^ [^mim ^max]}- 

Here, Ki are triangle elements. An example of such a mesh, built from the profiles of 
the different regions is shown on Fig. 9. Recall that the entire domain is meshed even 



inside the catch and collar regions because equations are set in the entire domain 
by the means of the permeability of the different media. 




Fig. 9. Unstructured mesh of the domain (10978 vertices - 21862 triangles) 

Mesh refinements are located near the region G, since it is the region where most of 
the turbulence occurs. 

The space discretization of the problem is based on the finite elements method. The 
velocity and pressure unknowns are approximated using P2/P1 finite elements. 

The associated discrete finite element spaces are the following: 



et Vr\TindK, = 0}, 



Qhi^r,z) = {qh e C\Qr,z),^K, G %,qH G P\K,)}. 



(103) 
(104) 



The turbulent kinetic energy k is approximated by P2 finite elements. 



The associated discrete finite element space is: 



m{nr,z) = {Wh e C°(fi,,,),Vir, G %,WH G P^K,), 



WhlvindK, = 0,Wh\r,ndK, = 0}. 



(105) 



The discrete weak formulation of the problem (31) — (36), (95) — (96) — (97) is the 
following: 



Finding {uh,Ph,h) e Vhi^r,z) x Qh{^r,z) x yVh{^r,z) tel que : 



[ -TT-^h \r\7rdrdz + I (u/iV)u/iV/i \r\'Kdrdz — f phV ■ Vh \r\Tidrdz 



+^ / ('^o + '^OlVu^ + (Vu;,)*) : (Vv;, + (Vv;,)*) \r\T^drdz 
+ / V{uh){t,:x.)v \r\TTdrdz 

- / V ■ uq\r\ndrdz = 0,yvh eVhi^r,z),yqh ^ Qhi^r,z); 
/ — — |r|7r(ir(i2; + / (u/iV)/i;/j |r|7r(ir(iz 
+ / i>t{Vkh : Vw/,) |r|7r(ir(i2; - / ^|Vu/, + (Vu/,)*pw/, |r|7rrfrrf2; 

+ --A-klwh\r\TTdrdz = Oywh eWhi^r,z)- 



(106) 



6.3 Time discretization 



Denote 6t the time step. Let u™, and k"^ be the time approximates of the mean 
velocity, the modified pressure and the turbulent kinetic energy respectively, at the 
time t"* = mSt. 



The convective terms in the problems are approximated using a characteristic Galerkin 



method [22], [13 



Consider a convective term like (uV)b. 

A Taylor expansion of the derivative 

Dh dh 



yields the approximation 



Dh ^ b"^+^ - (b'"(a; - u"'{x)St) 



(10^ 



Let X{x, t; s) be the solution of the problem: 

r dX 



ds 



u{X, s) 



(109) 



X\.=t = X. 



X{x,t; s) is the position at time s of the particule situated at position x at time t. 



Then: 



Dh b'^+i - b^oX" 



Dt 



6t 



(110) 



where X'"(a;) = X(x,r+^t"). 



Following 22|, an implicit scheme (see equation (111)) is chosen for the Navier- 
Stokes problem with eddy viscosity and a half-implicit one (see equation (112)) for 
the turbulent closure equation. 



For all m = 0, ... 



T 



5f 



find (ur\pr\C') e V,(fi.,,) X Q,(a,,.) X W,(^].,.) such as: 



^/ in'^+^ -u'^oX]^)wh\r\'Kdrdz- I p^^^^ V ■ Wh\r\'Kdrdz 
dt Jr^ JTh 

+\ I (z/o + Ci/v^)(Vu^+i + (Vu™+i)*) : (Vv^ + (Vvft,)*)|r|7rrfrrfz 

+ f V{vi';^^^)wMTxdrdz- ( V ■ Ml'+Uh\r\TTdrdz 

- I 0.0000001 \r\-Kdrdz = 0,Vv;, e V/,(a,.), Vg^, G Q/.(l).,.); 



[111) 



r 1 

6t Jr. 



{k^'+^ - k'f^oX]^) Wh \r\7cdrdz 



f {C2lVk^){Vk"'+^ : Vwh) \r\Tidrdz 



:il2l 



4 / {-CjVk^)\Vu^ + {Vu 



h.rn+1 
m\t\2 '^h 



Wh \rWdrdz 



+ 



'KK^^Wh\r\T^drdz = Q, 



where X™(x) is a numerical approximation of The parameters Ci (i=l, 2, 

3) are adimentionahzed constants. 



The penahzation term in equation (111) 

g,, 0.0000001 \r\'Kdrdz 



leads to a more regular problem [13 . 



The initial values for the velocity and pressure unknowns (u5J,p^) are obtained by 
solving an auxiliary Stokes problem, and the turbulent kinetic energy /cJJ is initialized 
to a constant in the entire domain. 



The solving process is iterative. As soon as the final time T is not reached, one 
solves numerically the kinetic energy problem, then the Navier-Stokes/Brinkman 
with eddy viscosity part, the time step is increased, the kinetic energy part is solved 
again and so on. 



6.J^ Parameters settings 



Different parameters have to be set to perform the simulations: 

- the parameter £(x) in the definition of the eddy viscosity function (see equation 
(15)) is defined as a constant in each triangle, its value in a triangle being equal to 
the longest edge of this triangle, 

- the water kinematic viscosity uq in equation (see equation (15)): uq = 1.141*10"^ 
mV^ at 15 °C, 



- the initial turbulent kinetic energy equal to a constant in the entire domain and 
equal to 0.01 m^s~^, 



- the permeability K in the different regions: 
10'' s^^ in the fluid region Vt^] 

10~^ s~^ in the catch region Gc and collar region G/; 
1 s~^ in the net region G^; (113) 

5 s~^ in the net region G^; 

6 s~^ in the net region G\. 

The unit of K is [s~^] since it is formally the ratio between the kinematic viscosity 
[m^s"'^] under a permeability surface [m^]. 

Simulations have shown that the subdomain G\ could be considered as permeable 
(i.e. as a fluid part). In fact, the mesh opening in G\ is so high that the meshes do 
not disturb the flow. 

- the time step set equal to 0.667 s, 

- the adimentionalized constants, found numerically: Ci = 0.1; C2 = 0.05; Ga = 
0.03. 



6.5 Numerical results 

Using the parameters defined in the previous section, we use the free software 
PreeFem++ to compute the fluid problem. Runs were made on a bi-processor Pen- 
tium Xeon EM64T 3.2Ghz, with 2Go RAM. 

The global behavior of the flow is shown on Fig. 10 where the streamlines are drawn. 

The level curves of the z component of the mean velocity are given in Fig. 11. 

Fig. 12 gives the level curves of u^. and the Fig. 13 and 14 gives those for the turbulent 
kinetic energy k. 

The use of an unstructured mesh leads to a slight asymmetry in the graphics for the 
turbulent kinetic energy. 



X(x) = 



Fig. 10. Streamlines 




Fig. 11. Level curves of behind the catch 
Those figures give several results: 

- a laminar fiow at the output (see Fig. 10). It allows us to keep the simplified 
boundary conditions (95) — (96) at the output. 



- the escapement of the inner velocity inside the net takes place just in front of 
the catch (see Fig. 12), 



- the turbulence is mainly located behind the catch and is low in the surroundings 
of the net (see Fig. 13 and 14), 



- two main eddies are located behind the catch (see Fig. 10, 11, 14). 




Fig. 13. Level curves for k in the surroundings of the net 

Let us compare now the experimental profiles given at the beginning (Fig. 3), mea- 
sured by a LDV technique for Uz with those obtained numerically (see Fig. 15-16-17). 




Fig. 14. Level curves for k behind the catch 
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Fig. 15. Profiles 2 and 3 after 50 iterations 




Fig. 16. Profiles 4 and 5 after 50 iterations 



One can see that the numerical profiles fit well with those obtained experimentally 
(see Fig. 15-16-17). 




An interesting feature is emphasized by computing the norm 2 of the difference of 
the velocity and the turbulent kinetic energy between two successive iterations (see 
Fig. 18). A stationary state is reached after about 50 iterations: the residual for u 
is equal to 0.00109346, and the one for k equal to 0.000406185. This is in agreement 
with the fact that we are studying mean quantities. 

To conclude, we have a model that leads to remarkable results in comparison with 
the available experimental data. In this particular case of a rigid net, our model 
looks appropriate. Moreover, this model has the advantage that its application to a 
3D problem is possible, especially if we make use of a fictitious domain technique 
that does not require a complex mesh generation. 
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